前言
这一部分是《Data Analysis for the life sciences》的第5章线性模型的第1小节。这一部分的主要内容涉及矩阵的表示方法。
背景知识
数据分析中的很多模型都可以用矩阵代数来表示。我们指的这些模型都是线性模型(linear models)。线性(Linear)并不是说它们是直线，而是说是线性组合。使用矩阵代数来描述线性模型非常方便，因为我们可以更好地构建这些模型，并使用数学工具来方便计算。在这一部分中，我们将详细地描述我们如何使用矩阵代数来构建并拟合线性模型。
在本书中，我们关注的线性模型是基于二分法的线性模型：例如治疗组与对照组。在前面提到的小鼠饲料的案例中就是这种线性模型的一个典型案例。在这一部分中，我们会描述一些比较复杂的案例，不过还会继续关注二分法变量的线性模型。
当我们学习线性模型时，需要记住，我们使用的变量仍然是随机变量。这就意味着，我们使用线性模型获得的估计值也是随机变量。虽然数学上这理解起来比较复杂，但是我们前面学到的一些概念在这里仍然适用。现在可以先做一些练习，在线性模型的练习中复习一些随机变量的概念。
练习
P172
设计矩阵(The Design Matrix)
R中的formual()和model.matrix()可以用于生成各种线性模型的设计矩阵(design matrices)（也被称为模型矩阵，即model matrices)。例如，在小鼠的饲料实验中，我们可以这样构建模型，如下所示：
其中，$Y_{i}$表示体重，而$x_{i}$等于1的时候，表示小鼠$i$吃的是高脂饲料(hf)。我们可以使用实验单位(Experimental unit)来表示N个不同的实验动物，这里的实验单位是小鼠。我们称这类变量为指示变量(indicator variables)，因此这类变量仅用于标明实验单位，表示某种干预的有无。现在用矩阵代数来描述就是以下形式：
或为：
简化形式为：
其中$\mathbf{X}$是设计矩阵，一旦我们定义了设计矩阵，那么我们就能计算出最小二乘估计。这个过程称为拟合模型(fitting the model)。在R中，我们可以直接在lm()函数中输入一个公式(formula)直接进行拟合。在后面的一个脚本中，我们会使用model.matrix()函数，这个函数是在lm()函数内部使用的。这种使用方式哦可以方便地让我们将R的formula与矩阵$\mathbf{X}$连接起来，有助于我们对lm()的结果进行解释。
设计方案
设计矩阵的选择是线性模型中的关键步骤，因为它编码了哪些系数将适合用于建模，以及样本之间的相互关系。一个常见的误解就是，选择实验设计要遵循简单明了的原则，即要在描述中直接标明样本。事实并非如此。关于每个样本（无论是对照组还是处理组，实验批次等）的基础信息里并不意味着这是一个“正确”的设计矩阵。设计矩阵还应该对$X$中的变量是如何解释$Y$的值这些信息进行编码，这是研究者们必须要考虑的。
在这一部分的案例中，我们使用线性模型来比较不同组之间的差异。因此，最终进行计算的设计矩阵应该至少包含两列，即截矩(intercept)列，截矩列是第1列，另外就是第2列，第2列指定样本。在这种情况下，线性模型中的两个系数(coefficient)分别为：第1个系数是截矩(interccept)，它代表了第1组的总体均数，第2个系数是差值，这个差值表示的是第2组与第1组总体的差值。当我们进行统计学分析时，往往第2个系数是我们最感兴趣的，因为我们想知道这两组数据是否有差异。
我们可以在R中通过2行代码来编制实验设计。我们先是使用波浪线~来构建一个公式。这种形式的公式意味着，我们想使用波浪线右侧的变量的观测值来构建统计模型。然后我们再放进一个变量的名称，就能知道样本来源于哪些组。
看一个案例。
现在我们有两组数据，分别是对照组和高脂饲料组(hf)小鼠的体重数据。现在我们用1来表示对照组小鼠，2来表示hf组。我们首先要告诉R，这些值不能被当成数值，而是应该当成因子。然后使用公式~group来构建模型，如下所示：
|
|
结果如下所示：
|
|
这里面有很多信息，attr后面的信息不用管。
其实上面的计算过程也可以使用model.matrix()来构建，如下所示：
|
|
如下所示：
|
|
如果我们不对组别(group)进行因子转换，那么就会出现下面的情况：
|
|
即：
|
|
这种形式不是设计矩阵。因为没有经过factor转换的数值就是单纯的数值，而非指示变量(indicator variable)，因子类型变量的数值与具体的数字无关，字符串也能生成这些变量，如下所示：
|
|
结果如下所示：
|
|
多组设计
现在我们设计了3组数据，如下所示：
结果如下所示：
|
|
现在我们的设计矩阵中有了第3列，这个第3列属于第3组。另外一种方法就是通过在公式中指定+0来添加第3组，如下所示：
结果如下所示：
|
|
注：在这一处我不太理解，原文是This group now fits a separate coefficient for each group.翻译为汉语则是这一组现在适用于每组的单独系数，等后文中有介绍的时候再补上，以下是其它的参考书（数学建模：基于R.薛毅 著）对这一内容的理解：
lm()中的参数formula是模型公式：
y ~ 1+x1 + x2这种形式，表示常数项，$X_{1}$，$X_{2}$的系数；y ~ x1 + x2这种形式，去掉了1，其意义与原来加1的一样。y ~ 0 + x1 + x2这种形式，添加了0，表示拟合成就齐次线性模型（或者是改为y ~ -1 + x1 + x2）。
多变量设计
前面我们构建的线性模型只有1个变量，即饲料。而在生命科学中，在一次实验中，往往会涉及多个变量，例如我们也许会研究饲料与性别这两个因素对体重的影响，在这种情况下，我们也许就会设计4组，如下所示：
组别如下所示：
|
|
假如我们认为饲料对体重的影响和性别对体重的影响（这只是假设）是相同的，那么我们就可以构建一个线性模型，如下所示：
现在我们在R中对这个模型进行拟合，此时我们需要添加另外一个变量来构建一个设计矩阵，如下所示：
如下所示：
|
|
这个设计矩阵包含一个截矩(intercept)，一个饲料变量(diet)，一个性别变量(sex)。我们可以说这个线性模型解释了组间差异与条件变量的差异。但是，如上所述，这个模型成立的前提是，我们假设了饮食和性别的差异会对小鼠的体重产生相同的影响。我们称这些情况为累加效应(additive effect)。对于每一个变量，我们添加一个效应时，不用考虑其它的变量是什么。针对这种情况，还有一种线性模型，这种线性模型里还会考虑两个变量之间的潜在交互作用。
如果我们要考虑交互作用的话，模型的构建如下所示：
|
|
结果如下所示：
|
|
指定比较水平relevel()
在构建设计矩阵时，我们所选水平的参考水平用于比较。默认情况下，我们都是按照第1个水平因素进行比较，但是我们可以使用relevel()函数来指定哪个水平因素作用为对照，例如我们选择第2组作为参考水平，如下所示：
|
|
结果如下所示：
|
|
这里涉及了relevel函数，此函数信息如下所示：
relevel()函数
全称reorder levels of factor，函数的功能是：对一个因子的水平重新排序，以通过指定ref来指定第1个水平作为参考水平（例如在回归中使用二元解释变量，可以指定哪个水平作为参考水平），其余的往后移。在contr.treatment对照中非常有用，contr.treatment通常会将第1个水平当作参考水平。
用法：relevel(x,ref,...)，其中x是一个未排序的因子，ref：指定的参考水平，通常是一个字符串。
或者是在factor()中直接明确地指出，如下所示：
|
|
结果如下所示：
|
|
model.matrix的计算
model.matrix()函数会捕获R全局环境中的变量，除非这些数据已经作为数据框的形式提供给到model.matrix()函数，如下所示：
|
|
结果如下所示：
|
|
连续型变量
前面内容中的变量都是指示数据（有的统计学书中指的是分类变量，哑变量），而在一些实验设计中，我们研究的是数值型变量（连续型变量），不需要将它们转换到第1列。例如在自由落体实验中，时间是一个连续型变量，如下所示：
|
|
结果如下所示：
|
|
在这个模型里，我们用到了I()函数，I()的功能是，对一个变量进行数学变换，在上面的公式里，I(tt^2)的意思是就是，通过I()，将tt^2这个变量当作是一个变量，如果不加I()，则是单纯地计算tt^2的值。查看I()的文档就可以发现，这个函数的全称是Inhibit Interpretation/Conversion of Objects。
在生命科学中，我们还有可能会涉及这样一种情况，例如我们想研究一个药物的量效关系，例如研究这个药物在0mg，10mg和20mg时的效果。
将连续型数据强制作为变量很难进行分析。这就是为什么指示变量只是假设2组的均值不同，而连续型变量会假设预测测变量和结果之间存在着某种特定的关系。
在自由落体实验中，我们有重力加速度这个理论的支持。而在父子身高案例中，由于这些数据是二元正态分布变量，如果我们加上一些条件的话，父子之间的身高存在着线性关系。但是，我们发现，如果不对年龄等这类变量进行“调整(adjust)”的话，并将它们包含在线性模型中，其统计结果就值得商榷。我们不支持将这样的变量加到线性模型中，除非已经有支持这种数据的模型出现。
练习
P183